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Spatially confined rigid membranes reorganize 
their morphology in response to the imposed con- 
straints. A crumpled elastic sheet presents a com- 
plex pattern of random folds focusing the defor- 
mation energy [Ij while compressing a membrane 
resting on a soft foundation creates a regular pat- 
tern of sinusoidal wrinkles with a broad distribu- 
tion of energy [2H8]. Here, we study the energy 
distribution for highly confined membranes and 
show the emergence of a new morphological insta- 
bility triggered by a period-doubling bifurcation. 
A periodic self-organized focalization of the defor- 
mation energy is observed provided an up-down 
symmetry breaking, induced by the intrinsic non- 
linearity of the elasticity equations, occurs. The 
physical model, exhibiting an analogy with para- 
metric resonance in nonlinear oscillator, is a new 
theoretical toolkit to understand the morphology 
of various confined systems, such as coated ma- 
terials or living tissues, e.g., wrinkled skin [3j, 
internal structure of lungs [10], internal elastica 
of an artery [llj, brain convolutions [12, 13j or 
formation of fingerprints [14j. Moreover, it opens 
the way to new kind of microfabrication design 
of multiperiodic or chaotic (aperiodic) surface to- 
pography via self-organization. 

Several theoretical approaches have been proposed to 
describe the wrinkling instability for very small compres- 
sion ratio, i.e. near the instability threshold [21 O [7]. 
However, the large compression domain remains largely 
unexplored with the notable exception of the wrinkle to 
fold transition observed in Ref. [8 for elastic membrane 
on liquid and the self-similar wrinkling patterns in skins 
[9]. In the former case, the deformation of the membrane 
is progressively focalized into a single fold, concentrating 
all the bending energy. In contrast, for thin rigid mem- 
branes on elastomers, large compression induces pertur- 
bations of the initial wrinkles but the elasticity of the soft 
foundation maintains a regular periodic pattern whose 
complexity increases with the compression ratio. 

A PDMS film, stretched and then cured with 
UV/ozone, or a thin polymer film bound to an elastomer 
foundation, remains initially flat. Under a slight com- 
pression, S = (Lo - L)/Lo, these systems instantaneously 
forms regular (sinusoidal) wrinkles with a well-defined 
wavelength, Aq. Increasing S generates a continuous in- 
crease of the amplitude of the wrinkles and a continuous 



shift to lower wavelength (A = Xo(l - 5) see Fig. [T^). 
By further compression of the sheet, more complex pat- 
terns emerge. Above some threshold, S > S2 - 0.2, we 
observe a dramatic change in the morphology leading 
to a pitchfork bifurcation: one wrinkle grows in am- 
plitude at the expense of its neighbours (Fig. [T]). The 
profile of the membrane is no longer described by a sin- 
gle cosinusoid but requires a combination of two periodic 
functions, cos and cos The amplitude of the 2A 
mode increases with the compression ratio, while the A 
mode vanishes. This effect is similar to period-doubling 
bifurcations in dynamical systems jT5] |16] observed in, 
for example, Rayleigh-Bernard convections [17], dynam- 
ics of the heart tissue [T8H2Q] , oscillated granular mat- 
ter 122 ^ or bouncing droplets on soap film [23]. In 
contrast to previous works, we describe here a spatial 
period- doubling instability which is rarely observed [24] . 
Nonlinear coupling between two modes, one with double 
the wavelength of the other, also appears in post-buckling 
of cylindrical shells as reported in the classical work of 
Koiter (see [25 and references therein). 

The thin inextensible membrane of length Lq is com- 
pressed horizontally by a distance A = Lq - L along the 
X-axis and is bound to an elastic foundation that ini- 
tially fills the half-space ^ < 0. The system is assumed 
to remain invariant in the z direction (see Fig. [T]). The 
projected length along the x-axis, Lq - A, is given by 

Lo-A= r^°d^cos(^, (1) 
Jo 

where £ is the arc length measured along the curve. The 
quantity (j) is the angle between the tangent to the surface 
and the horizontal. The derivative of this angle with 
respect to the arc length, 9^^, gives the local curvature 
of the membrane (for clarity, partial derivatives such as 
^ are written as di). The relative compression ratio is 
given hy 6 = A/Lq. 

The response of this thin membrane resting on an elas- 
tomer substrate is determined through minimization of 
the energy per unit of width, U. Two energetic contribu- 
tions are to be considered: i) the elastic bending energy 
of the thin sheet, 

ub = ^ l"^" dm<k)\ (2) 

where the parameter Bm is the bending stiffness of the 
membrane {Bm Emh^ ^ Em being its Young's modu- 
lus and h its thickness); ii) the energy of deformation 
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FIG. 1: Evolution of morphologies, wavelengths and amplitudes with compression. System 1 (a,b): A PDMS 
foundation is cured with UV/ozone which modifies the elastic properties of its surface. The thickness of the membrane is about 
10 - 20 /xm depending on the irradiation time. The wavelength, Ao, of the initial wrinkling instability is about 50 - 100 /xm. 
System 2 (c,d,e): A thin coloured stiff PDMS film is bound to a thick soft PDMS foundation. The thickness of the membrane 
is about 200 fim and the initial wrinkle wavelength is about 3 mm. The compression ratios, 6, are equal to 0.165, 0.19 and 0.24 
for panels c, d and e respectively, (f) The systems are compressed uniaxially along the x-axis. The wavelength and amplitudes 
of the wrinkles are measured for successive values of the relative compression 6. (g) Amplitudes (Ai, A2) and wavelength. A, 
as a function of the compression ratio 6. Experimental data for system 1 are reported with symbols •, ■ and ♦ for 30 min., Ih 
and 2h of irradiation respectively whereas the symbol v is used for system 2. Results of the linear (dotted lines) and nonlinear 
(solid lines) theories are also reported. Before period-doublin g, th e expression of the amplitude A computed from Eqs. ([T]) is: 
^ = ^ (1 - f - w)- The wavelength A is computed from |ll|: ^ = 1 - ^ + 0(^^ (B/XofS^), where B is the amplitude of 
the subharmonic mode. 



of the elastomer. The constraint of inextensibility of the 
membrane ([T]) is taken into account with the help of a 
Lagrangian multiplier F identified to the cross-sectional 
pressure per unit length. The Euler-Lagrange equation 
obtained from the energy of the system gives the equilib- 
rium of normal forces along the membrane and is given 

by 

B^dty + Fdly + Py^O, (3) 

where y and Py are functions describing the vertical el- 
evation of the membrane and the normal pressure from 
the elastomer acting on the membrane, respectively. At 



linear order, Py = Kl-L{diy)^ where K is the stiffness 
coefficient of the foundation proportional to its Young's 
modulus {K = 2E{1 - cr)/(l + cr)(3 - 4cr), where a is the 
Poisson ratio, see Supplementary Information) and H is 
the Hilbert transform. The first nonlinear contribution 
due to the elastomer can be computed for periodic defor- 
mation with one mode of frequency q and Eq. ([3| reduces 
then to 

Bmdjy + Fd^.y + Kqy + K2q^y^ = 0, (4) 

where K2 is also proportional to the Young's modulus 
{K2 = £;(l-2cr)(13-16a)/2(l + cr)(3-4cr)^ see Supple- 



3 



mentary Information). Notice that, as for the hnear re- 
sponse of the susbtrate, the nonhnear term involves also 
an Hilbert transform for multimode profile. Due to the 
quadratic nonlinearity from the foundation, the equation 
[4] giving the profile of the membrane implies an up-down 
symmetry breaking: vertical extension and compression 
along the y-axis are no longer equivalent. This equation 
can also be viewed as a spatial equivalent of a nonlinear 
oscillator, like a simple pendulum, with which it shares 
many similarities. 

Equation Q reduces to a linear oscillator for small 
amplitudes of the instability. In this regime, the period 
is independent of the amplitude, as for a simple pendu- 
lum, in agreement with observation and usual theories. 
Indeed, nonlinear terms can be neglected for small am- 
plitudes and the curvilinear and Cartesian coordinates 
coincide: ^ x, (j) ^ dxV- Equation Q admits sinusoidal 
solutions y{x) = Acos(27rx/A) provided that the pressure 
F and the wavelength of the wrinkling instability are re- 
lated by 



(5) 



This relation shows that below a threshold F < Fc = 
SqQBm there is no associated wavelength and the mem- 
brane stays flat. At the threshold, F = Fc, the wrin- 
kling instability emerges and a unique and constant wave- 
length, Ao, is selected ^3^ i26j 



(6) 



The selection of this particular wavelength is obtained 
from a minimization of the energy through a minimiza- 
tion of F. The inextensibility constraint ([T]) gives the 
evolution of the amplitude of the instability as a func- 
tion of the relative compression, A = ±Ao\/^/7r. However, 
neither the evolution of the wavelength with 6 nor the 
period-doubling bifurcation are captured by this linear 
model. 

To determine the supercritical morphology, we study 
the stability of the single wavelength pattern in the 
weakly nonlinear regime. We thus consider a small pe- 
riodic perturbation, eu, characterized by a frequency /c, 
of the nonlinear solution for the shape of the membrane: 
y ^ y + eu, e being arbitrarily small. The equation for the 
perturbation, u, in the leading order in the amplitude A 
of the instability is then given by 



Bmdju + Fd'^u + Kku ■■ 



-2K2kAqo cos{qoi)u. (7) 



The term appearing in the right-hand side of this equa- 
tion is due to the quadratic nonlinearity of the founda- 
tion (stiffness K2). Interestingly, this equation is similar 
to the Mathieu equation, describing resonance in para- 
metric oscillators [27H3Q] . For usual forced oscillators, 
like a simple pendulum with a variable length (the most 
famous example of this resonance is given by the giant 



censer, O Botafumeiro [29]), the unforced system is char- 
acterized by a given period and the additional frequency 
needed to produce a parametric resonance is provided by 
an external agent. For all amplitudes of the forcing, the 
resonance appears provided that forcing and oscillator 
frequencies are related through k = qo/2. 

In our system however, we should also consider a con- 
straint related to the minimization of F {i.e., minimiza- 
tion of energy since U{S) = Lq /q F{5')d5' where 5 is 
the relative compression) determining the amplitude of 
the forcing term at which the 2A mode emerges. Actu- 
ally, the period-doubling instability cannot be observed 
for amplitudes smaller than a threshold {i.e., defining a 
compression threshold, (^2). 

From equation we can deduce that the profile 
should be described by a multimode solution of the form, 
y{^) - Efcli C'/c cos(/c(7o^/2). Indeed, without any loss of 
generality, the wrinkled pattern can be assumed to be 
described by an even function since the system is invari- 
ant under horizontal translation. The numerical analy- 
sis of equation Q, adapted to multimode periodic solu- 
tions, shows a very good agreement with experimental 
data (Fig. [l^). Notice that the convergence is already 
reached with the four first modes (see Supplementary In- 
formation) . The relevance of the model is further demon- 
strated by the excellent agreement between experimen- 
tal and calculated profiles (Figs. ^ and d). We should 
emphasize that the model relies on a single parameter, 
K2IK, that determines the period-doubling threshold 82- 

In order to preserve an explicit analysis and to cap- 
ture the physics of the model, we restrict the following 
discussion to the ansatz y{t) = Acos{qoi) + B cos{qoi/2). 
Substituting this ansatz in Eq. Q, we obtain a system 
of two equations in A, B and F, admitting two solu- 
tions. A trivial solution corresponds to the evolution 
before period-doubling: F^q^Bm) = 3 and B = {A be- 
ing determined by the inextensibility constraint). The 
second solution involving a subharmonic mode (5^0) 
reads 



F 
B' 



17/4 + 2A, 
2A(5 + 8A). 



(8) 
(9) 



where F = F/(q^Bm), A = 2K2Aqo/K and B = 
2K2Bqo/K. Equation (39) is no longer invariant under 



a change of sign of A. Indeed, the amplitude A of the 
harmonic mode can be either positive or negative since 
the nonlinear system is characterized by a up-down sym- 
metry breaking due to the quadratic nonlinearity of the 
foundation upon deformation. Fig. [2^ and b shows the 
evolutions of both solutions with the amplitude, A. The 
symmetry breaking induces two regimes. For A> 0, F is 
always larger than the value associated to the harmonic 
mode alone, i.e. F = 3. The corresponding shape for 
the membrane is forbidden and thus not observed ex- 
perimentally, see Fig. (2)3. In contrast, for A < 0, the 
emergence of a subharmonic mode is energetically fa- 
vorable {F < 3) beyond a threshold value, see Fig. [2^. 
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FIG. 2: Predictions of the model and comparison with experimental profiles, (a), (b) Evolution of the cross-sectional 
pressure F and of the amplitude B of the subharmonic mode as a function of the amplitude A of the harmonic mode. Arrows 
indicate the path followed by the systems during compression. (Insets) Representative shapes of the membrane according to 
the value of F. (c) Comparison between theoretical and experimental profiles of the membrane for system 1 for two values of 
5. (d) Comparison between theoretical and experimental profiles of the membrane for system 2 for two values of 5. 



From Eq. ([9|, we observe that B starts to grow precisely 
from this threshold. This analysis does not, however, 
imply that an harmonic mode with a positive ampli- 
tude, A > 0, is stable against subharmonic perturbations. 
Indeed, the above analysis is performed using, without 
loss of generality, an even function to describe the evo- 
lution of the wrinkled pattern. Having found the ener- 
getically favorable pattern in this case, we can use the 
translation invariance to generate equivalent patterns: 
y{l - tt/qo) = -Acos{qoi) + Bsm{qoi/2). The sign of A 
being now reversed, it implies that an harmonic mode 
with a positive amplitude is also unstable against sub- 
harmonic perturbations above the same thershold and 
leads to the same wrinkled pattern but translated. 

Through the inextensibility constraint, the threshold 
for A implies the existence of a critical relative compres- 
sion, (^2, for the onset of the period-doubling instability. 
Using the relation between A and S at the lowest order. 



we obtain 

The critical compression needed to observe a period- 
doubling bifurcation for wrinkling instability, 62^ strongly 
decreases with the Poisson ratio of the elastic foundation. 
The values found numerically for the ratio K2/K ~ 0.25 
yield a Poisson ratio around 0.44 which is close to the 
value usually reported in literature for PDMS (~ 0.48). 

Moreover, this model based on nonlinear oscillator 
should imply that, for larger amplitudes of the 2A mode, a 
period-quadrupling bifurcation characterized by a wave- 
length 4A would appear. This behaviour is indeed ob- 
served in Fig. |3^, b for compression ratios larger than 
0.26. This last observation clearly suggests that cascades 
of spatial period- doubling bifurcations can be observed 
for the elastic instability of rigid membrane, provided 
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FIG. 3: Additional consequences of the up-down symmetry breaking on wrinkled patterns, (a) Evolution of ampli- 
tudes for large compressions showing the period-quadrupling bifurcation in system 2. The solid curve is obtained numerically 
for K2IK = 0.27. The dashed curve is added to help visualizing the second bifurcation characterizing the period-quadrupling 
instability, (b) Wrinkled structure showing the period-quadrupling instability {S ^ 0.37). (c) Profile of a thin stiff PDMS 
membrane resting in between two identical soft PDMS foundations for a relative compression 6 ^ 0.23 (d) Profile of a thin stiff 
PDMS membrane bound to a soft PDMS foundation for a relative compression ^ 0.23. 



that the up-down symmetry is broken. Such a cascade is 
known to lead to chaos after several bifurcations [15f Il6j. 
There is however a geometric limitation in our system in 
contrast to previously reported temporal period-doubling 
cascade. Indeed, the evolution of the pattern saturates 
as soon as sharp folds appear (see Fig.jsfD). For instance, 
due to finite thickness of the membrane, we experimen- 
tally reached at most period-quadrupling structures. 

A further confirmation of our approach can be ob- 
tained. Our interpretation of the period-doubling bi- 
furcation in rigid membrane on elastomer implies that 
the dynamics should be governed by nonlinear terms of 
even order, which break the up-down symmetry Conse- 
quently, systems with an up-down symmetry, like a thin 
elastic membrane resting on a liquid [8] , do not develop a 
period-doubling instability. Interestingly, we could make 
trilayers restoring the symmetry. Indeed, a system com- 
posed of a thin elastic membrane in-between two identical 
soft foundations, one below and one above the membrane, 
does not exhibit the period-doubling bifurcation. Instead 
it develops patterns similar to those observed with float- 
ing membranes. In Fig.jSfD and c, we compare the profile 
of the membrane when there are one or two foundations 
for the same compression ratio. 

The second salient feature of the nonlinear wrinkling 
instability is the continuous decrease of the wavelength 
with the compression ratio 5. This effect arises from the 
change from curvilinear to Cartesian coordinate. The 
wavelength is measured along the horizontal x-axis while 
the shape of the membrane is determined in curvilinear 
coordinates £ where it is constant. For a periodic profile 
y{i), with a wavelength A^, A^^ = A is given by 



The evolution of the wavelength along the horizontal x- 
axis at the leading order in the amplitude of the insta- 
bility, A, is given by 



A 1 (All! 

Ao " A2 



1-S 



(12) 



A= / dicos(l)= / diJl 
Jo Jo 



(11) 



in very good agreement with experimental data in 
Fig. [5. 

The universal model describing the formation of 
wrinkled patterns based on nonlinear oscillator dynam- 
ics should explain observations in very different fields. 
For example, a better understanding of the elastic 
instability of rigid membranes will help to determine 
the exact mechanisms leading to the growth of wrinkled 
morphology in living systems. It is also a new blueprint 
to develop multiple-length-scale microfabrication tech- 
niques useful in the design of specific topography. 

Methods 

Experiments were carried out using polydimethylsilox- 
ane (PDMS) elastomer (Sylgard 184) purchased from 
Dow Corning. Two different systems were studied. Sys- 
tem 1: a bare elastomer of PDMS is irradiated with 
UV in presence of oxygen. Ozone is generated and will 
affect the crosslinks density of the PDMS outer sur- 
face. The rigidity of the surface drastically increases 
with the irradiation time to finally yield a brittle over- 
layer covalently bound to the uncured elastomer. Sys- 
tem 2: multilayers prepared by a simple assembly of 
monolayers of different elastic properties. The "rigid" 
and "soft" layers correspond to elastic modulus values of 
1200 and 10 kPa, respectively. To ensure a very strong 
adhesion between both PDMS films and avoids delam- 
ination during the compression, these two PDMS elas- 
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tomers were assembled by contact after a plasma cur- 
ing (in a Plasma Cleaner oven). The experimental set- 
up was a custom-built stretching/compressing device. 
The UV/03 modified PDMS was compressed by using 
a stretched/curing/release experiments. The measure- 
ments were achieved using image analysis from micro- 



tomed slices of the samples. The bilayer PDMS assem- 
bly were compressed by inducing a macroscopic radius 
of curvature. The measurements were performed from 
macro photography of the cross-section of the samples 
(see Supplementary Information for further details). 



Witten, T. A. Stress focusing in elastic sheets. Rev. Mod. 
Phys. 79, 643-675 (2007). 

Bowden, N., Brittain, S., Evans, A. G., Hutchinson, J. 
W. & G. M. Whitesides. Spontaneous formation of or- 
dered structures in thin films of metals supported on an 
elastomeric polymer. Nature 393, 146-149 (1998). 
Cerda, E. & Mahadevan, L. Geometry and physics of 
wrinkling. Phys. Rev. Lett. 90, 074302 (2003). 
Vandeparre, H. et al. Slippery or Sticky Boundary Condi- 
tions: Control of Wrinkling in Metal-Capped Thin Poly- 
mer Films by Selective Adhesion to Substrates. Phys. 
Rev. Lett. 99, 188302 (2007). 

Vandeparre, H. & Damman, P. Wrinkling of Stimulore- 
sponsive Surfaces: Mechanical Instability Coupled to Dif- 
fusion. Phys. Rev. Lett. 101, 124301 (2008). 
Huang, J. et al. Capillary wrinkling of floating thin poly- 
mer films. Science 317, 650-653 (2007). 
Jiang, H. et al. Finite deformation mechanics in buck- 
led thin films on compliant supports. PNAS 104, 15607- 
15612 (2007). 

Pocivavsek, L. et al. Stress and fold localization in thin 
elastic membranes. Science 320, 912-916 (2008). 
Efimenko, K. et al. Nested self-similar wrinkling patterns 
in skins. Nature Mater. 4, 293-297 (2005). 
Diamant, H., Witten, T. A., Ege, C, Gopal, A. & Lee, 
K. Y. C. Topography and instability of monolayers near 
domain boundaries. Phys. Rev. E 63, 061602 (2001). 
Strupler, M. et al. Second harmonic microscopy to quan- 
tify renal interstitial brosis and arterial remodeling. J. 
Biomed. Opt. 13, 054041 (2008) 

Richman, D. P., Stewart, R. M., Hutchinson, J. W. & 
Caviness, V. S., Jr. Mechanical model of brain convolu- 
tional development. Science, 189, 18-21 (1975). 
Toro, R. & Burnod, Y. A Morphogenetic Model for the 
Development of Cortical Convolutions. Cereb. Cortex 15, 
1900-1913 (2005). 

Kiicken, M. & Newell, A. C. A model for fingerprint for- 
mation. Europhys. Lett. 68, 141-146 (2004). 
Feigenbaum, M. J. Quantitative universality for a class 
of nonlinear transformations. J. Stat. Phys. 19, 25-52 

(1978) . 

Feigenbaum, M. J. The universal metric properties of 
nonlinear transformations. J. Stat. Phys. 21, 669-706 

(1979) . 

Libchaber, A., Laroche, C. & Fauve, S. Period dou- 
bling cascade in mercury, a quantitative measurement. 
J. Physique 43, L211-L216 (1982). 

Guevara, M. R., Glass, L. & Shrier, A. Phase locking, 
period-doubling bifurcations, and irregular dynamics in 
periodically stimulated cardiac cells. Science 214, 1350- 
1353 (1981). 

Fox, J. J., Bodenschatz, E. & Gilmour, R. F. Period- 
doubling instability and memory in cardiac tissue. Phys. 



Rev. Lett. 89, 138101 (2002). 
[20] Berger, C. M. et al. Period-doubling bifurcation to al- 
ternans in paced-cardiac tissue: Crossover from smooth 
to border-collision characteristics. Phys. Rev. Lett. 99 
058101 (2007). 

[21] Melo, F., Umbanhowar, P. B. k Swinney, H. L. 
Hexagons, kinks, and disorder in oscillated granular lay- 
ers. Phys. Rev. Lett. 75, 3838-3841 (1995). 
[22] Venkataramani, S. C. & Ott, E. Spatiotemporal bifurca- 
tion phenomena with temporal period doubling: patterns 
in vibrated sand. Phys. Rev. Lett. 80, 3495-3498 (1998). 
[23] Gilet, T. & Bush, J. Chaotic bouncing of a droplet on a 

soap film. Phys. Rev. Lett. 102, 014501 (2009). 
[24] Losert, W., Shi, B. Q. & Cummins, H. Z. Spatial period- 
doubling instability of dendritic arrays in directional so- 
lidification. Phys. Rev. Lett. 77, 889-891 (1996). 
[25] Hutchinson, J. W. & Koiter, W. T. Postbuckfing theory 

Appl. Mech. Rev. 23, 1353-1366 (1970). 
[26] Groenewold, J. Wrinkling of plates coupled with soft elas- 
tic media. Physica A 298, 32-45 (2001). 
[27] McLachlan, N. W. Theory and application of Mathieu 

functions, Dover, New- York, 1962. 
[28] Blanch, G., Chapter 20: Mathieu Functions, in Handbook 
of Mathematical Functions with Formulas, Graphs, and 
Mathematical Tables, Eds. Abramowitz, M. & Stegun, I. 
A., Dover, New- York, 1972. 
[29] Sanmartm, J. R. O Botafumeiro: Parametric pumping 

in the Middle Ages. Am. J. Phys. 52, 937-945 (1984). 
[30] Van den Broeck, C. & Bena I. Parametric Resonance 
Revisited, in "Stochastic Processes in Physics, Chem- 
istry, and Biology", Eds. Freund, J. A. & Poschel T. Lect. 
Notes Phys. 557, 257-267 (2000). 
Acknowledgements The authors thank T. Witten, 
B. Davidovitch, H. Diamant, S. Desprez, C. Troet- 
sler, S. Gabriele and G. Carbone for fruitful discus- 
sions. This work was supported by the Belgian Na- 
tional Funds for Scientific Research (Mandat Impul- 
sion Scientifique), the Government of the Region of 
Wallonia (CORRONET and REMANOS Research Pro- 
grammes) and the European Science Foundation (Eu- 
rocores FAN AS programme, EBIOADI collaborative re- 
search project). F.B. acknowledges financial support 
from a return grant delivered by the Federal Scientific 
Politics. 



7 



Supplementary information for "Multiple-length- 
scale elastic instability from period-doubling bi- 
furcation cascade mimics parametric resonance of 
nonlinear oscillators" 



I. MATERIALS AND METHODS 

Experiments were carried out using polydimethylsilox- 
ane (PDMS) elastomer (Sylgard 184) purchased from 
Dow Corning. By changing the proportion of crosshnker, 
we have adjusted the elastic properties of the elastomer. 
The elastic moduli measured for various crosslinked 
PDMS elastomers is given in Fig. [4] [1]. 

Two different systems were studied: i) UV/03 cured 
PDMS films (system 1) and, ii) multilayers prepared by 
a simple assembly of monolayers of different elastic prop- 
erties (system 2). 

By irradiating a bare elastomer of PDMS with UV 
in presence of oxygen, we generate ozone molecules that 
modify the crosslinks density of the PDMS outer surface. 
As demonstrated by numerous studies [2], the rigidity 
drastically increases with the irradiation time to finally 
yield a brittle overlayer covalently bound to the uncured 
elastomer. As showed below, the thickness of modified 
PDMS layer increases with irradiation time. 

In contrast, the multilayers systems were prepared 
by assembling two PDMS films of different rigidity 
(high/low crosslinks density with high/low elastic modu- 
lus). The "rigid" and "soft" layers correspond to elastic 
modulus values of 1200 and 10 kPa, respectively. After a 
plasma curing during four minutes (in a Plasma Cleaner 
oven), these two PDMS elastomers were assembled by 
contact. The plasma curing ensures a very strong adhe- 
sion between both PDMS films and avoids delamination 
during the compression of the multilayer. 

The experimental set-up was a home-made stretch- 
ing/compressing device. The UV/03 modified PDMS 
was compressed by using a stretched/curing/release ex- 
periments (maximum compression ratio, 5 ^ 0.55). The 
measurements were achieved from image analysis from 
microtomed slices of the samples perpendicular to the 
symmetry axis of the pattern. Due to the very small 
wavelength observed for UV/03 systems (10- 100 /im) 
a microscope was required. The bilayer PDMS assem- 
bly were compressed by inducing a macroscopic radius 
of curvature (maximum compression ratio, S ^ 0.4). The 
measurements were performed from macro photography 
of the cross-section of the samples. 



II. INFLUENCE OF UV/03 ON PDMS 
ELASTOMERS 



As demonstrated by several studies, UV/03 curing of 
PDMS generates a rigid overlayer by a chemical trans- 
formation of the PDMS chains, the chemical composition 
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FIG. 4: Evolution of the elastic modulus of the PDMS elas- 
tomers as a function of the fraction of crosslinker (the fraction 
recommended by Dow Corning is 10% to obtain the more rigid 
elastomers) . 
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FIG. 5: Evolution of the wavelength at the threshold, Ao, with 
the curing time in UV/03 oven (irradiation time). The solid 
line corresponds to a fit of the data with a rational scaling 
law (Ao 



for very long irradiation time corresponding to brittle sil- 
ica [2]. In addition, the wrinkling instability can be used 
to check the finite size of the modified PDMS surface 
(Fig.|§. 

Considering the relation Aq oc h{Eral E)^!^ (where /i, 
Em and E are the thickness, the elastic modulus of the 
rigid "membrane" and the elastic modulus of the founda- 
tion, respectively), the evolution of the initial wavelength 
with irradiation time is in agreement with a diffusive be- 
haviour. Indeed, the UV/03 curing affects the outer sur- 
face of the PDMS film up to an effective thickness, /leff, 
expected to grow with the irradiation time according a 
\/t law. 



III. PARAMETRIZATION, EQUATION FOR 
THE MEMBRANE AND ORIGIN OF THE 
NONLINEARITY 

The membrane is assumed to be inextensible and in- 
variant along one direction such that it is completely de- 
termined by its {x^y) Cartesian coordinates. Let ^ be 
the curvilinear coordinate and (j) the angle between the 
tangent to the surface and the horizontal direction. Ac- 
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cording to elementary differential geometry, the shape of 
the surface is then parametrized by 



i^^y) = I j ds cos (I), I ds sin 



, J ds sin ^ j , 

= (13) 

where / gives the height of the membrane as a function 
of the curvilinear coordinate. In addition, the derivative 
di(j) is equal to the curvature of the membrane. 

The thin rigid membrane resting on the soft thick foun- 
dation is assumed to be invariant along one direction. 
Consequently, the evolution of the system under com- 
pression is described in a plane (x^y). For small de- 
formations, the equation giving the shape of a free and 
inextensible plate resting initially ai y = and invariant 
along one direction reads [3] 



y 



{IV) 



(14) 



where B^n - £^m^^/(12(l - cr^)) is the flexural rigidity. 
Em is the Young modulus, a is the Poisson coefficient and 
h is the thickness of the plate. F is the in-plane applied 
load on the plate. For large deformations, the nonlinear 
terms are of odd order and the first nonlinearity is cubic. 
Indeed, this system has an up-down symmetry: bending 
up or down a plate is energetically equivalent. 

If an elastomer fills the semi-plane ^ < 0, an additional 
term, describing the force normal to the surface due to 
the elastomer, must be added to Eq. (14). For small 



deformations, the linear theory of elasticity can be used. 
The equation satisfied by the displacement vector, is 
given by [3] 



(l-2cr)Ai2 + V(v-i2) =0, 



(15) 



where cr is the Poisson ratio. We assume that there 
is no deformation along the z-axis; the deformation is 
then studied in the (x, ?/)-plane. Under this condition, 
Eq. (IT5| reduces to 



2{l-(j)dlu, + {l-2(j)dlu, + d: 
2{l-(j)dluy + {l-2(j)dlu ' 



xy^y 



3 V 



= 0, 
= 0. 



(16) 



The boundary conditions are Ui{x^ y -oo) = (z = x, ?/), 
Ux{x^y = 0) = and Uy{x^y = 0) = /(x), where f{x) is 
a known function describing the shape of the surface of 
the foundation. The vertical force per unit area at the 
surface of the elastic substrate subject to deformation 
f{x) is given by [3] 



Py{X,y = 0) = (JyiUi 



(17) 



where there is summation on repeated indices and where 
ft is the unit normal vector to the surface. At the lowest 
order in the amplitude of /(x), we simply have n = (0, 1) 
and Py is given by (Jyy{x^ y = 0). To compute this quan- 
tity, we solve Eq. ([16]) with the appropriate boundary 



conditions given above. Once the displacement vector is 
known, the strain tensor is computed, in linear theory of 
elasticity, from [3] 



1 / dui duk \ 
Uik = - - — + — — . 

2 \ oxh oxi I 



The stress tensor is then given by 

E I cr 



l + cr 



(18) 



(19) 



To solve Eq. ( 16 ), we use the Fourier transform of Ui{x^ y) 
along the x-axis 

Ui{k,y) ^ T{ui{x,y)) ^ —= / Ui(x,y)e"'''dx. (20) 



The PDF system ( 16 ) reduces then to this ODE system 

(1 - 2a) dyUx - ikdyUy - 2(1 - a) k'^Ux = 0, 
2{l-(j)dluy-ikdyUx-{l-2(j)k^Uy = 0. (21) 

This last system of equations can be decoupled to obtain 

0, 



dyUx - 2k^dlux + k'^Ux 



ikdyUy ■ 



2{l-(7)k^Ux-{l-2(7)dyUx = 0. (22) 



2 ~ 



The solution which satisfies the boundary conditions 
reads 



'^y{k,y) 



(3 -4a) 

m 

(4<7-3) 



(4a - 3 + |fc|t/) el'^l^ (23) 



Using Eqs. (18) and (19) in Fourier space, we obtain 



Py{k,y = 0) = K\k\f{k), 



where 



K 



2{1 -a)E 
(l + cr)(3-4cr)' 



(24) 



(25) 



E being the Young's modulus. Consequently, the quan- 
tity to be added to Eq. ( 14 ) is given by 



Pyix,y = 0) = KJ^-\\k\fik)). 



(26) 



This expression still involves the Fourier transform of 
f{x). It is possible to obtain an equivalent form involving 
only f{x). The Hilbert transform, of a function f{x) 
is a linear operator given, by definition, by the convolu- 
tion of f{x) and h{x) = l/{7rx): n{f{x)) = {h f){x). 
Using the properties of the Fourier transform we then 
have 



TinidJix))) = ^2T:T{h)T{dJ{x)), 



= sgn(k)kfik) = \k\fik). 



(27) 
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Using this last relation, we can now write 



(28) 



Due to the property of the Hilbert transform (hke 
l-L{cos{qx)) = sm{qx) and 1-L{sm{qx)) = -cos(gx)), for 
a periodic profile characterized by only one frequency 
this additional term is equivalent to qy. Consequently,in 
the linear regime, the equation to solve for a periodic de- 
formation characterized by only one frequency q is then 
given by 



Br, 



(IV) 



+ Fy'' + Kqy = 0, 



(29) 



where K = K{E^ a) is the stiffness of the elastic founda- 
tion. 

Now, we extend this result for the response of the sub- 
strate to the case of larger deformations where nonlin- 
earities are not negligible. Here, we compute the first 
nonlinear correction to the response of the foundation 
restricting our analysis to deformations relevant to our 
system. The purpose of this paper is to understand the 
emergence of the subharmonic mode. We thus perform 
a weakly nonlinear analysis valid up to the threshold of 
the period-doubling instability. 

For larger deformations, the contribution of the elastic 
foundation must be computed using nonlinear elasticity 
theory. The relation between the strain tensor Uik and 
the displacement vector u is now given by [3] 



W^ + ^ + ^^'j (30) 
2 \ dxk dxi dxi dxk J ' 



where there is summation on repeated index. In linear 
theory of elasticity, the quadratic term on the right-hand 
side of (30) is neglected (see Eq. ([l8|). It is however 
possible to rederive the equation of elasticity taking into 
account this nonlinear term since the equation satisfied 
by the displacement vector is obtained at equilibrium by 
daik/dxk = 0, where the expression of the stress tensor is 
given by (19). This nonlinear equation reads 



(1 

+ (l-2a)^ 



d'^Ui dui 



u 



(31) 



M 



dxl dxi 



^ d'^U£ dui 
dxkdxi dxk 



0, 



with i = 1,2,3. For a given shape of the deformation of 
the elastic foundation, the normal force per unit area, P, 
is still computed from the relation Pi = aikUk, where n is 
the unit normal vector to the surface. Once P is known, 
it can be added to (29). Setting rj = \/5, all quantities are 



expanded in power of r]. Equation (|31| is then solved in 
perturbation with appropriate boundary condition and 
in particular Uy{x^y = 0) = Aco^{qx) + C cos{2qx) ^ where 
A is assumed to be of order 77 and C of order rf . The 
nonlinear terms of Eq. ( [31] ) being quadratic, the first non- 
linear mode of the wrinkled pattern of our system will 
be of the form cos{2qx). This nonlinear shape used as 
boundary condition should describe accurately enough 



the profile of our system up to the emergence of the sub- 
harmonic mode. Sufficiently near the threshold of the 
period-doubling instability, the amplitude of this sub- 
harmonic mode remains very small with respect to the 
dominant mode. The nonlinearity associated to this sub- 
harmonic mode should thus remain negligible. 

The lowest nonlinear response of the elastic defor- 
mation to a periodic deformation of the system with 
frequency q is then found to be Py = Kq{Acos{qx) + 
2Ccos(2(7x))+ {K2/2)A'^q^cos{2qx), with K2 = E{1 - 
2cr)(13 - 16cr)/2(l + a)(3 - 4a)^ At order r/^ this ex- 
pression can be written as Py = Kqy + K2q^{y^ - (^^}), 
where (•} = f^-dx. For our purpose, since we are 
considering periodic solutions, it can also be written as 
Py = Kn{y') + K2VHW? - iW?)]' To simplify the 
notations, the average term is dropped in the following 
equations, and in the main text, but is taken into account 
in the computations. The resulting nonlinear equation, 
valid only for wrinkled pattern characterized by only a 
single mode of frequency q^ is then found to be 



S^t/^'"^' + Fy" + Kqy + K2qV = 0. 



(32) 



In this equation, cubic or higher order nonlinear terms, 
coming from the bending of the membrane or the defor- 
mation of the foundation, are neglected. Notice that for 
multimode profile, characterized by several frequencies 
qi, the nonlinear equation involving the Hilbert operator 
should be used. 



IV. PERIOD-DOUBLING BIFURCATION 

We consider the equation that the shape of the surface 
of the system, must satisfy at the quadratic order 



+ Fy" + Kn{y') + K^Hiy'f = 0. 



(33) 



The solution at quadratic order that minimizes the en- 
ergy, i.e. minimizing F, can be written as 

y = Acos{qo£) + super-harmonic terms of order 77^, (34) 

with A = ±277/^0, given by the inextensibility constraint 
and T] = \/S. We perform now a linear stability analysis of 
this solution against periodic perturbations, character- 
ized by a single frequency k. We substitute the function 
y-\-eu into Eq. ( 33 ) , e being an infinitesimal quantity. At 



the first order in e and r] we obtain 

BmU^^^^ +Fu" + Kn{u') = -2K2n{y')n{u), (35) 

= -2K2AqQCOs{qQ£)n{u'). 

Since u is characterized by a single frequency /c, this equa- 
tion can be written as 

BmU^^^^ + Fu' + Kku = -2K2Aqok cos{qoi)u. (36) 

This equation is identical to Eq. (7) reported in the 
main text. We rescale now the parameters, the indepen- 
dent variable and the function as follow: F = Fq^Bm^ 
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K2 = K2qQBm, Qo^ = z and = u/qq. Equa- 



K 

tion ( |36[ ) takes now the form 

ii^^^^ + F^x'' + 2a;ii = -2AqoK20J cos{z)u, (37) 
where uj = k/qo. Assuming u = cos(cjz), we obtain 
(uj^ - Fuj + 2)cos(ujz) = -AqoK2[cos((l-uj)z) 



+ cos((l . 



(38) 



If uj < 1/2, keeping the lowest order Fourier modes, the 
equation is satisfied provided 2/uj. However for 

these values of cj, F is always larger than the value ob- 
tained for the harmonic mode, i.e. F = 3. Consequently, 
this mode, cos(cjz), cannot emerge because the shape 
adopted by the system is the one that minimizes F. If 
UJ > 1/2, the equation cannot be satisfied by keeping only 
the lowest order Fourier mode. However, if u = 1/2, the 
equation is satisfied if 



F = 17/4 + 2AqoK2. 



(39) 



The amplitude of the harmonic mode, A, is either pos- 
itive or negative since the system, at the linear order, 
has an up-down symmetry. However, adding a quadratic 
nonlinearity breaks this symmetry, this is why Eq. (39) 
is no longer invariant under a change of sign of A. For 
large enough (in absolute value) negative values of the 
amplitude A, i.e. -AqoK2 > 5/8, Eq. (39) shows that 



adding a subharmonic mode to the shape of the mem- 
brane leads to a smaller value of F than the value ob- 
tained with an harmonic mode alone. This critical am- 
plitude implies a threshold, S2, for the period-doubling 
instability through the inextensibility constraint. This 
mechanism not only explains the emergence of the sub- 
harmonic mode but also selects the correct sign for A 
leading to profiles actually observed in experiments, see 
Fig. [6j Obviously, this analysis does not, however, imply 
that an harmonic mode with a positive amplitude, A > 0, 
is stable against subharmonic perturbations. Indeed, the 
above analysis is performed using, without loss of gen- 
erality, an even function to describe the evolution of the 
wrinkled pattern {y{i) = Acos{qo£)-\- B cos{qo£/2)). Hav- 
ing found the energetically favorable pattern in this case, 
we can use the translation invariance to generate equiv- 
alent patterns: y{i - ir/qo) = -Acos{qoi) + B sm{qoi/2). 
The sign of A being now reversed, it implies that an 
harmonic mode with a positive amplitude is also unsta- 
ble against subharmonic perturbations above the same 
threshold and leads to the same wrinkled pattern but 
translated. However, the amplitude of the subharmonic 
mode cannot be obtained with the linear equation (37) 
and is computed in Sec. [Vl 



V. NUMERICAL ANALYSIS OF THE MAIN 
NONLINEAR EQUATION 

In order to confirm the period-doubling mechanism 
proposed in the main text, we analyze numerically the 
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FIG. 6: Schematic views of the possible shapes of the mem- 
brane after period-doubling according to the sign of the har- 
monic mode. 



nonlinear equation (5) in the main text, namely 

B^y^'""^ + Fy" + Kn(y') + K2Hiy'f = 0. (40) 



We use a rescaling similar to the one proposed in Sec. IV 
except for the function: F = Fq^Bm^ K = 2qQBm^ 
K2 - K2qQBm^ qo^ = ^ and y = u/{qoK2). Equation 



(40) reduces then to 

^(^^) + Fu'' + 2n{u') + n{uf = 0. 
We search solutions under the form 

N 

u{z) = ^ c/c cos(/cz/2), 

k=l 



(41) 



(42) 



with Ck = 2qo{K2/K)Ck where Ck are the expansion co- 
efficients introduced in the main text. Substituting the 
form (|42| into ( [41] ) we obtain a system of N equations 
with + 1 unknowns {F and q, i = 1, . . . ,N). N un- 
knowns can be expressed in terms of C2, for example. 
This last coefficient is determined from the inextensibil- 
ity constraint. 

For N = 2^ the computation is analytical and is already 
sufficient to understand the mechanism for the emergence 
of the subharmonic mode. Remembering that the terms 
u and u'^ are operators that should be multiplied by the 
frequency of the modes on which they act, we obtain the 
following system of equations 



16 



C2(3-F) + ^ = 
(17-4F + 8C2) = 0. 



(43) 



This system admits two solutions. The first one is F = 3 
and Ci = (c2 being determined by the inextensibil- 
ity constraint). This solution correspond to the evolu- 
tion of the shape of the membrane without subharmonic 
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FIG. 7: Evolution of F and the amplitudes A 1,2 (see Fig. 1 of 
the main text for the definition) as a function of the rescaled 
amplitude of the harmonic mode for several values of N. 



mode. The second solution characterized by a subhar- 
monic mode is 



17 



+ 2C2 



cj = 8c2(F-3) = 2c2(5 + 8c2). 



(44) 
(45) 



The expression for F is identical to Eq. (|39|) obtained 
in Sec. 



IV 



The amplitude of the harmonic mode, C2, 
can be either positive or negative since the system, at 
the linear order, has an up-down symmetry. However, 
adding a quadratic nonlinear ity breaks this symmetry, 
consequently the expression (44) for F is no longer in- 



variant under a change of sign of C2. For large enough 
(in absolute value) negative values of the amplitude C2, 
i.e. -C2 > 5/8, Eq. (44) shows that adding a subhar- 
monic mode leads to a smaller value of F than the value 
obtained with the harmonic mode alone, i.e. F = 3. The 
situation is summarized in Figs. 2a and b in the main 
text. 

For = 4, only the solution without subharmonic 
mode can still be obtained analytically in a simple form; 
it reads 



F = 4-yi + 

ci = cs=0 
2c4 = 1-A/1 + 



m 



(I) 



(46) 



We see again that the quadratic nonlinearity breaks the 
up-down symmetry since for any sign of C2, C4 is always 
negative. The solution containing a subharmonic mode 



C: 1e- 




FIG. 8: Evolution of the Fourier mode coefficients as a func- 
tion of the rescaled amplitude of the harmonic mode for N = 8. 



or solutions for larger N are obtained numerically. In 
Fig. [t) we study the convergence for the evolution of F 
and Ai^2 (see Fig. 1 of the main text for the definition) 
when the solution is expanded as Eq. (42) for several val- 



ues of A^. Convergence is essentially reached for = 4. In 
Fig. [8) we present the evolution of the coefficients of the 
Fourier modes as a function of the rescaled amplitude, 
C2, of the harmonic mode for N = 8. We find that the 
subharmonic mode emerges for 02=02- -0.42. Return- 
ing to the original variable, we find that the subharmonic 
mode emerges when the amplitude, A, of the harmonic 
mode reaches the value 



A 



0.42 



4n(K2/K) 



(47) 
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FIG. 9: Threshold, S2, of the period-doubling instability as a 
function of the strength of the nonlinearity K2/K. The low- 
order relation (48) (blue dashed curved) is compared to an 
higher order relation obtained from Eqs. (47) and (49) (black 
solid curve). 

The amplitude of the harmonic mode is related to the 
relative compression 6 through the inextensibility con- 
straint. We thus obtain an expression for the threshold 
in compression, 62^ for the emergence of the subharmonic 
mode. At the lowest order, we have A/Xq = ^/S/tt which 
leads to 



0.105 



XK2/K)) 



(48) 
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in Fig. [9| 

Finally, in Fig. [Toj we present the numerical evolution 
of the pressure F as a function of the relative com- 
pression 5 together with corresponding profiles of the 
membrane. 



FIG. 10: Numerical evolution of F as a function of the rel- 



ative compression 5. The analytical expression (46) for the 



evolution of F without subharmonic mode is also presented. 



An higher order expression A/Aq as a function of 5 reads 

(49) 



Ao TT I 8 



128 



This last expression together with (47) yield an higher 



order relation for 82 which is plotted, together with (48), 
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